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We define a new set of excitations in the XY model which we call "fractional vortices". In the 
frustrated XY model containing -n bonds, we make the ansatz that the ground state configurations 
can be characterized by pairs of oppositely charged fractional vortices. For a chain of tt bonds, 
the ground state energy and the phase configurations calculated on the basis of this ansatz agree 
well with the results from direct numerical simulations. Finally, we discuss the possible connection 
of these results to some recent experiments by Kirtley et al [Phys. Rev. B 51, R12057 (1995)] on 
high-T c superconductors where fractional flux trapping was observed along certain grain boundaries. 
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I. INTRODUCTION 

The classical XY Hamiltonian is one of the most stud- 
ied models in statistical physics. In its usual, unfrus- 
trated form, it is written 



H 



(ij) 



1 



i((f>i ~ <t>j)], 



(1) 



where 4>i is a phase variable on the i th site (0 < < 2w) , 
the sum runs over distinct pairs (ij), and Jy is the en- 
ergy of the coupling between sites i and j. In the fer- 
romagnetic, nearest-neighbor case, Jjj vanishes except 
between nearest-neighbor sites and all the 's are equal 
to a single positive constant J. In this case, for spatial 
dimensionality d > 3, there is a phase transition to a 
ferromagnetic state at a critical temperature, with con- 
ventional critical phenomena. If d = 2, there is instead 
the Kosterlitz-Thouless-Berezinskii phase transition, in 
which pairs of oppositely charged integer vortices unbind 
at a finite temperature Tktb&- The classical XY model 
has been found to describe a wide variety of systems with 
complex scalar order parameters, including bulk super- 
conductors in d = 3, superconducting films, Josephson 
junction arrays in d = 2, and superfiuid He4 filmsB. 

Recently, the XY model with anti ferromagnetic bonds 
i.e, with bond strengths Jijrfk^ (also called it bonds), 
has received much attentionETel, in particular due to its 
possible relevance to hkrli^T c superconductors and other 
experimental systemst2rt3. Specifically, if we consider 
the grain boundary between two high-T c superconduc- 
tors with suitable misorientation of the crystalline axes, 
then the resulting Josephson coupling across the bound- 
ary can have the coupling energy Jy < 00. This is a 
consequence of the d x 2_ y 2 symmetry of the order param- 
eter in many high-T c materials. Such grain-boundary 
interfaces have lately been studied in a variety of ex- 
periments and in several geometries. These experiments 
have led to interesting results, such as the observation of 



the trapping of half-integer and also other fractional flux 
quantaliJ~tL3. These resujte_r;an be explained using mod- 
els involving it bondsuE3~Llj. Similar models involving it 
bonds have also been developed to explain such phenom- 
ena as the paramagnetic Meissner effectO, also observed 
in samples of high-T c superconductors. 

A key concept in understanding the effects of it bonds 
is "frustration." Consider, for example, the XY model on 
a square lattice with only the nearest-neighbor couplings 
nonvanishing. If a plaquette has an odd number of bonds, 
that plaquette is frustrated, in the sense that no choice 
of angles in the four grains making up the plaquette can 
simultaneously minimize all the bond energies. Thus, a 
single 7T bond will cause the two plaquettes adjoining that 
7r bond to become frustrated. In a line of tt bonds, only 
the two plaquettes at the end of the line will become 
frustrated. Because of the frustrated plaquettes, it is 
non-trivial to find the ground state of the XF-model with 
7r bonds. In this paper, we will show, both numerically 
and by analytical arguments, that these ground states 
are characterized by certain spatial phase configurations 
which we call fractional vortices. We will also derive an 
expression for the interaction energy of two fractional 
vortices in the XY model. 

The rest of the paper is organized as follows. In Section 
II, we define the fractional vortices and calculate the in- 
teraction energy of a bound pair of fractional vortices for 
the XY model. In Section III, we study the ground state 
of XY-lattices containing a single tt bond, two tt bonds, 
and a string of tt bonds. In each case, using a varia- 
tional ansatz for the ground-state configuration, we find 
that there is a critical tt bond strength, above which the 
ground state contains pairs of oppositely charged frac- 
tional vortices. To check these results, we directly cal- 
culate the ground-state energy of these lattices using a 
numerical relaxation technique based on the equations 
of motion for overdamped Josephson junctions. We find 
that both the ground-state energy and the critical it bond 



1 



strength, predicted by the variational approach, are in 
excellent agreement with the numerical results. Finally, 
in Section IV, we discuss the possible relevance of these 
numerical results to experiments carried out in systems 
containing tt junctions, such as high-T c superconductors 
containing grain boundaries, and tricrystals, as recently 
studied by Kirtley et aO~H 



II. FRACTIONAL VORTICES IN THE 
UNFRUSTRATED XY MODEL 

Consider the Hamiltonian (|l|) for an XY model defined 
on a square lattice with N x N sites. If all the nearest- 
neighbor couplings are equal, this may be written 



H= J^[l-cos( 

(ij) 



(2) 



Hereafter, we shall use units such that J = 1. The phase 
angle, at point (xi,yi) due to a fractional vortex of 
charge q at point (xq, yo) is defined to be 



(t>i{xa,yo,q) = q x tan 1 ( 



Vi-yo 



x 



Xj 



(3) 



For q = 1, we recover the standard configuration for an 
integer vortex. This definition can be seen as a gen- 
eralization of the concept of half-vortices introduced by 
Villained for the same model. Note that, while for the 
integer vortex the bond angles change continuously, the 
fractional vortex case is characterized by a branch cut, 
across which the bond angles are discontinuous. 

This singularity leads to several other distinctions be- 
tween integer and fractional vortices. For example, the 
energy associated with a single integer vortex is propor- 
tional to hi(N). In the thermodynamic limit, this is a 
weak divergence which makes the KTB vortex-antivortex 
unbinding transition possible. By contrast, the energy of 
an unbound fractional vortex is oc N, since the number of 
bonds along the branch cut is oc N. Thus, it is energet- 
ically unfavorable at all temperatures to create isolated 
fractional vortices. But a bound pair of fractional vor- 
tices with charges q and — q is much less expensive ener- 
getically, because then the branch cut is restricted to the 
line joining the two charges: the total energy should be 
proportional to the separation of the fractional vortices. 
For fixed q and large enough separations, this energy is 
always larger than that of a pair of oppositely charged 
integer vortices, whose energy varies as the logarithm of 
their separation. Nevertheless, for fixed separation, it is 
always possible to find a non-integer q such that that the 
energy of the fractional vortex pair is less than the corre- 
sponding energy for the integer vortices. In the following, 
we derive expressions for the energy of a bound pair of 
fractional vortices in the XY model, and compare them 
to numerical results obtained by calculating the energy 
explicitly for these configurations. 
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FIG. 1. Calculated energy of a bound pair of integer 
vortices as a function of separation n, in units of J. 
Full curve: numerically exact results. Dot-dashed curve: 
KT-approximation [Eq. (j^)]. Asterisks: KT-approximation 
plus core correction. 

We first consider a bound pair of integer vortices of 
charge ±1 located at (xq, yo) and {x\, yi). The standard 
KT expression for the energy of the pair is obtained by 
approximating the Hamiltonian as 



H 



<ij> 



if 



(4) 



For the phase configuration, we use (f>i = 4>i{xo, yo, +1) + 
<pi(xi,yi,-l), where xi - x = n and y\ - y = Q (in 
units of the lattice constant a) . Substituting this config- 
uration into (^) gives the Kosterlitz-Thouless formula for 
the interaction energy, Ekt, of two oppositely charged 
integer vortices: 



E KT (n) = 2tt 



7I " 

inn H — 



(5) 



In Fig. [l], we compare this expression to the energy of a 
pair of oppositely charged integer vortices, computed us- 
ing the same phase configuration but the exact H. The 
discrepancy arises from the expansion of the cosine fac- 
tor, which is inaccurate for the bonds closest to the vor- 
tices. This inaccuracy is remedied by a making a core 
correction i.e, by calculating the contribution from the 
bonds on the perimeter of the plaquettes surrounding 
the vortices exactly, rather than by a quadratic expan- 
sion. For large n, the core-correction energy, E c {n), is 
approximately given by 



E c {n) 



12 



(6) 



2n 2 16n 4 ' 

As can be seen from Fig. |l|, the numerically calculated 
energy is well approximated by Exrin) + E c {n). 
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FIG. 2. Schematic drawing of a bound pair of fractional 
vortices arranged parallel to the x-axis. For this configura- 
tion, the vortices are separated by three plaquettes. <f>i U and 
<f>U are the phases at the two ends of bonds in the shaded 
region (region A); the angle Qi,A in the text is defined by 
9i,A = 4>iu — <t>u- Core corrections are calculated only for the 
bonds denoted C\ and C2, as discussed in the text. 

Next, we calculate the energy E(q, n) of a pair of frac- 
tional vortices ±3, separated by a distance n in the x- 
direction, as shown schematically in Fig. |2[ Note that 
E(l,n) is simply the interaction energy of a pair of in- 
teger vortices, as just discussed. The phase configura- 
tion of this pair is then given by 4n — <f>i(xo,yo>o) + 
4>i{xi, yi, — q), where x\ — xq — n and y\ — yo = 0. We 
now divide the bonds into two groups: (A) those inter- 
sected by the line segment joining the two vortex cen- 
ters; and (B) the remaining bonds. Let E A (q,n) and 
Eb (q, n) be the corresponding energy contributions to 
E(q, n) coming from these two groups of bonds. Thus 
E(q,n) = E A (q,n) + EB(q,n). To obtain E(q,n) we 
proceed as follows: 

(1) . We calculate E A (1, n), using the quadratic expan- 
sion for the cosine. E A (l,n) is simply the contribution 
to the total energy of a bound pair of integer vortices 
arising from the bonds along the branch cut. Note that 
there are n bonds in region A. Once E A (l,n) is known, 
we get E B (l,n) = E KT (n) - E A (l,n). 

(2) . We obtain EB(q,n) by noting that, in the 
quadratic approximation, Es(q,n) = q 2 Es(l,n). 

(3) . Finally, we determine E A (q,n) by directly evalu- 
ating it using the full expression for the cosine, not the 
quadratic expansion. This is necessary, because the bond 
angles in region A are not small for arbitrary q. 

We now use the outlined procedure to obtain E(q, n). 

Step 1: Let &i A {l,n) = <fii U — <f>u denote the i th bond 
angle (cf. Fig. ^) in region A for q = 1. For the two- 
vortex configuration, 0i iA (l,n) is given by 



i, A (l,n) = 2 



tan 



1 



2i-l 



-tan 



1 



2n-2i+l 



(7) 



Using the quadratic approximation, the corresponding 
energy contribution E A (l,n) is given by 



E A {l,n) = -Y.OlAihn). 



(8) 



Using the approximation: tan 1 [l/(2i — 1)] ~ 1 / (2z — 1) 
for i > 2, we find 



E A {l,n) = - 



2 2n- 1 



[7+21n2-2-2n] 



E B {l,n) = E KT {n) - E A {l,n), 



(9) 



where ip{ x ) an( i V^O*-) are t ne Digamma function and its 
derivative, 7 is Euler's constant, and Ekt^) is given by 
Eq. (|). 

Step 2. Using the results of step 1 and Eq. (H) , we get 

E B (q,n) =q 2 [2tt In n + tt 2 - E A (l, n)} . (10) 

Step 3. The next step is to calculate E A (q,n). Bond 
angles in region A are given by 

6 itA (q,n) = q[2w-9 i , A {l,n)], (11) 

Correspondingly, the energy E A (q,n) is given by 

n 

E A {q,n) =5^(1- cos [0 iiA (?,n)]). (12) 

i=l 

Since the bond angle 8i iA (q,n) is not small for an ar- 
bitrary q, we cannot expand the cosine term only to 
second order. But for any q, the difference 9i tA (q,n) — 
9n/2. A { < lT n ) is a small parameter for any i > 2. Expand- 
ing the cosine term in Eq. (|l2|) to second order in this 
parameter, we obtain an expression for E A (q,n). This 
expression can be summed, and eventually gives 



E A (q,n) = (n-2) 



( 1 — cos a n ) H 5- cos a n H sm a n 



2 <^ 1-cos 



n-l 



37T 

~~2 



2n-l 



4q 2 



n-l 



cos a n + q sm a n 



E m , A (hn) + ^-cosa n £ ^ >A (l,n), (13) 



m=2 



m=2 



where 

n-l 



5^0 m ,A(l,n) = 27-4 + 41n2 + 2^(n-i) , (14) 



m=2 
n-l 



EO(l,")=7r 2 + -(7 + 21n2-2-2n) 



m=2 



a n = q[ 2?r 



(15) 
(16) 
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FIG. 3. Energy of a bound pair of fractional vortices in an unfrustrated XY-lattice as obtained numerically (*) and from the 
analytical approximation, Eq. ([L7|) (full line), for (a) fixed charge, q = 0.8, as a function of separation; and (b) fixed separation, 
n = 50, as a function of charge. 



Adding up the contribution from the two regions, we 
finally get the required expression for the energy of a 
bound pair of fractional vortices. As noted earlier, the 
core corrections must be included to attain high numeri- 
cal accuracy. In the present case, it is sufficient to include 
these corrections only for the bonds labeled Ci and C2 in 
Fig. ^|, using the procedure outlined earlier. This approx- 
imation is equivalent to extending region A to include 
bonds Ci and C2 . Correspondingly, the core-corrected 
total energy is given by 



E(q,n)=E A (q, n)+E B (q, n) - q l Q l c + 2 [1 - cos(q0 c )] 



where 



2 



(17) 



(18) 



2n + 1 

Expressions ( |l7j ) and (|l^) are compared to the results 
of numerical computation in Figs. ||(a) and|^(b); agree- 
ment between the two is excellent. On the basis of this 
agreement, which is equally good for all values of q and 
n which we have considered, we present this result as 
a good analytical expression for the interaction energy 
between two fractional vortices in the unfrustrated XY 
model on a square lattice. This result is a generalization 
of the integer vortex excitations proposed by Kosterlitz 
and Thouless. 

For large n, we can further simplify the above expres- 
sion by dropping terms of 0(1/ n) and smaller in Eq. fll?| ) 
to get 

E(q, n) = (n- 2) [1 - cos(27r<7)] + 2 In n [nq 2 - q sin(27rg)] 



2„2 



-tt q 



1 — cos( ) 



(19) 



III. FRACTIONAL VORTICES IN THE XY 
MODEL WITH tt BONDS 

The fractional vortex configurations introduced in the 
previous section provide a natural way of characterizing 
the ground state of the XY model containing tt bonds. In 
this section, we implement this description by making a 
variational guess for the ground-state configuration using 
fractional vortices. We then compare our variational re- 
sults with those obtained by numerically relaxing to the 
ground-state configuration, and find excellent agreement. 
In the following subsections, we will focus on obtaining 
the critical bond strength, A c , above which the ferromag- 
netic ground-state solution becomes unstable, and the 
ground-state configuration contains bound pairs of frac- 
tional vortices. For the case of one and two tt bonds, we 
also compare our results to those from previous studies 
by Vannimenus et au. Note that in these calculations, in 
which the goal is to calculate the threshold bond strength 
above which the ferromagnetic ground state becomes un- 
stable rather than the absolute energies as a function of A, 
it is unnecessary to include the core corrections. Hence, 
A c can be calculated analytically as demonstrated below. 
The ground-state configuration and energy for A > A c 
do need the core corrections for greatest accuracy. We 
obtain them numerically using our variational guess and 
discuss them in the subsequent section. 



A. One tt bond 

Wc first consider the case of a single tt bond i.e, a single 
antifcrromagnetic bond in a host of ferromagnetic bonds. 
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FIG. 4. Schematic plot of three configurations, each containing -k bonds (full lines), which are (a) parallel and adjacent, (b) 
parallel and non-adjacent, and (c) perpendicular and non-adjacent. Also shown are the corresponding locations of the fractional 
vortex charges for the variational configurations. 



As before, we take the bond strength of the ferromagnetic 
bonds to equal unity, and we denoted the strength of the 
antiferromagnetic bond by A (A > 0). The problem is 
to obtain the ground-state configuration and energy for 
arbitrary strength, A, of the ir bond. 

To solve this problem, we make a variational guess for 
the ground-state configuration: it is the phase configura- 
tion corresponding to a bound pair of fractional vortices 
of strength ±q, located at the centers of the two plaque- 
ttes adjacent to the ir bond. The charge, q, is a varia- 
tional parameter with respect to which the ground-state 
energy is minimized for a given A. 

The total energy of this configuration discussed above 
is readily obtained using the procedure of the previous 
section, suitably corrected for the fact that we have a it 
bond instead of a ferromagnetic bond. The angle differ- 
ence across the ir bond is given by 9^ — qir. Then, using 
Eqs. (§) and @, we get 



E B (q) 



1 



2 2 

q ir , 



while from Eq. (fL2J) we find 



E A (q) = 1 + Acos(<2tt). 



(20) 



(21) 



Adding these two terms gives the total energy of the 
configuration. Minimizing this energy with respect to 
q yields the condition 



qTT = A sin((?7r). 



(22) 



For A < 1, the ground-state configuration corresponds 
to q = 0: all the phases are perfectly aligned. For A > 1, 



the ground-state configuration corresponds to a bound 
pair of fractional vortices with charges ±q obtained by 
solving Eq. (p2|). Thus, the ferromagnetic ground state 
is unstable above a critical bond-strength value A c = 1. 
The same value has been obtained previously by workers 
using different approachesQ'Q. 



B. Two n bonds 

We now consider the case of two parallel, adjacent ir 
bonds. As before, our variational guess for the ground 
state is the configuration corresponding to a bound pair 
of fractional vortices; we take these to be located as 
shown in Fig. f|(a). The corresponding total energy is 
again calculated using the procedure outlined in Sec- 
tion II. Using Eqs. (||) and (|l0|), the energy contribution 
E B (q) is 



E B {q)=q 2 I 2?rln2 



-2 tan" 1 (1/3) 



(23) 



Using Eq. (|l2|), we get 



E A (q) = 2 + 2Acos{2q [3tt/4- tan -1 (1/3)]} . (24) 

Adding these two contributions gives the total energy, 
which is to be minimized with respect to q for a given A. 
This procedure gives the critical value A c = 0.563, which 
is in good agreement with the-jexact value A c = n/2 — 1, 
obtained by Vannimenus et aa. 

Next, we consider the case of two parallel, but non- 
adjacent 7T bonds, as shown in Fig. 0(b). Taking the 
bond centers to have the coordinates (0,0) and (m,n), 
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FIG. 5. Schematic plot of the assumed variational configu- 
ration containing m pairs of fractional charges, for a chain of 
7r bonds (solid line segments) of length 2m. 

we calculate the energy using the variational procedure 
described above with vortex cha rges as sh own. For large 
separation between the bonds (\/m 2 + n 2 3> 1), this pro- 
cedure gives 



E B(q) = q 2 [27T 2 - 2a mn — (n — a mn ) 2 ] 



and 



E A (q) = 2 + 2A cos [<?(7r 



where 



(to 2 



,2)2 ■ 



(25) 



(26) 



(27) 



Minimizing the total energy gives the critical bond 
strength as 



1 - 2a mn /TT 
1 + 2a mn /n 



(28) 



Similarly, for two non-adjacent perpendicular tt bonds 
[Fig. |(c)], we find 



E B (q) = q 2 [2tt 2 - 2/3„ m - (tt - (i mn f 



and 



E A (q) = 2 + 2A cos [q(n + f3, n 



where 



ftn 



2mn 



(to 2 



In this case, the critical bond strength is 

1 - Pmn/n 



1 + /9mnA' 



(29) 



(30) 



(31) 



(32) 



These results are identical to those obtained previ- 
ously by Vannimenus et al, using a different approachu. 
The agreement lends support to our hypothesis that the 
ground-state configuration of such systems can be charac- 
terized by a set of fractional vortices (in the cases consid- 
ered here, a set of only two oppositely charged fractional 
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FIG. 6. Critical bond strength A c as function of chain 
length, for a chain of n bonds of length n. Open circles and 
full curve: numerical results. Asterisks and dashed curve: 
analytical approximation, A c = 

vortices). Besides having the merit of simplicity, our ap- 
proach also easily yields the ground-state configuration 
and energy for arbitrary A. Moreover, our procedure can 
be used to obtain the ground-state even when the sep- 
aration between the bonds is not large. In particular, 
for two parallel bonds such that to = n — 1, our varia- 
tional procedure yields the surprising result that A c — 1 
for this configuration. The same result was obtained in 
a numerical study done by Gawiec et aB. Finally, our 
variational ansatz can readily be generalized to longer tt 
bond chains, as we shall see in the next section. 



C. Chains of n bonds 

Next, we consider chains of tt bonds of length n > 3. In 
this case, we make the variational ansatz that the ground 
state consists of n/2 or (n + l)/2 pairs of oppositely 
charged fractional vortices for even or odd n, arranged 
as shown in Fig. [|. As before, we proceed by calculating 
the contribution to the total energy from regions A and 
B. However, the procedure outlined in Section II has to 
be generalized to include many pairs of fractional vor- 
tices. Since the details are significantly different from 
that outlined in section II, we briefly describe the gener- 
alized procedure below. 

(1). We consider the case in which all charges have 
magnitude unity. The total energy of this configuration 
is given simply by the KT expression 



E KT = -2tt 



j<k 



3 



(33) 
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FIG. 7. Ratio of i bond angle, 6i, to the first bond 
angle 9\, plotted as a function of i for chain lengths 
n = 5(A), 15(D), 25(*),35(o). The bond strength is A = 1. 
Inset: variation of 0i with chain length n, for A = 1. 

where njk is the distance between the charges q.j and qu 1 
and the second sum runs over all the individual charges, 
each of which has magnitude unity. This result is ob- 
tained by using a small-angle expansion for the contri- 
bution from each bond angle difference 9b, and summing 
those contributions to give 



Ekt — 



Id 



\2 
■>) 



The bond angle 9b is, in turn, decomposed as 
9b = ^ qkOk 



'k.b, 



(34) 



(35) 



where k labels the position of the charges and 9k,b is the 
contribution to 9b from a charge of unit magnitude at k. 

(2). Next, we consider the case in which the charges 
are fractional (\q\ < 1). The bonds can still be divided 
into classes A and B as discussed earlier. In the case of 
fractional charges, the bond angle differences in region B 
are still given by Eq. (|35|). Correspondingly, the energy 
contribution from bonds in region B is 



Er 



1 B 

5E' 



, A+B 

- y < 

2 V 

b 



(36) 



where ^ b and ^ b denote sums over all bonds in region 
B and in region A. 



FIG. 8. The i bond angle, 6i, plotted as a function of i 
for bond strengths A = 0.5, 1.0, 1.5,2.0. The chain length, n, 
is taken to be 20. 

(3). Finally we calculate the energy Ea, using suitable 
expressions for the bond angles in region A, as obtained 
from the multi-vortex configuration but without making 
the small angle expansion. 

To minimize the resulting total energy, which is a func- 
tion of all the qk's, we used two procedures: (i) Powell's 
multidimensional direction set methodE3, and (ii) a ge- 
netic algorithms. Both methods successfully converged 
to the same minimum energy and configuration, from 
which we deduced the critical bond strength, A c , for var- 
ious values of n. Fig. || shows our results for A c (n). As 
can be seen from the main part of the figure, it fits very 
well to the approximate expression A c ~ 1.16/n. A con- 
sequence of this 1/n dependence is: if the system has 
a finite concentration of 7r bonds randomly distributed 
in the lattice, then in the thermodynamic limit A c — * 0. 
This behavior follows from the fact that, in the ther- 
modynamic limit, there is always a finite probability of 
having an arbitrarily large chain length n, and hence an 
arbitrarily small A c . 

We now look at the variations in the bond angles along 
the 7r bond chain as a function of chain length, n, and 
bond strength, A (for A > A c ). First, we discuss the vari- 
ation with fixed bond strength, taking A = 1. Fig. |?] 
shows the ratio of the bond angles, 9i, along the chain 
(not including the central bond) to the bond angle across 
the first 7r bond, 6\, as a function of position along the 
chain for various chain lengths. A number of features 
deserve mention. First, for a chain of length n = 2m, the 
ratio of the bond angles, 6>;/6>i, for i < m is independent 
of m. Second, since the bond angle distribution is sym- 
metric, we only need to look at the bonds in the range 
1 < i < to. Third, the bond angles increase monotoni- 
cally as one moves along the chain from its edges to 
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FIG. 9. Bond angle, 8, as function of bond strength A for a 
single 7r bond obtained using (a) fractional-charge variational 
ansatz (*), and (b) numerical simulations (solid line). 

the center i.e, as i increases from 1 to m. The inset shows 
the variation of 0\ with chain length for A = 1. Note that 
with this choice of A, the bond angles saturate quickly: 
they are roughly constant over the "interior bonds," such 
that i > 3. Furthermore, this constant value (approxi- 
mated by the central bond angle) approaches tt as the 
chain length n increases. 

Fig. |H shows how the bond angles, Qi, vary with bond 
strength, A, for a fixed chain length (n = 20). As already 
seen in the previous figure, 6i rapidly tends to saturate 
towards its central value with increasing i. Moreover, 
the central bond angle quickly increases from to tt as 
A increases for A > A c . Thus, we can 'tunc' the central 
bond angle to any desired fraction of 7r by appropriately 
adjusting A. 

Although the underlying variables are the 9i's, it is of 
interest to mention corresponding trends in the fractional 
charges. For small A, these charges decrease monotoni- 
cally with increasing i, so that the largest charges reside 
at the ends of the chain. For larger A, charges compara- 
ble in magnitude to those at the ends appear away from 
the ends. 



D. Numerical Check of Variational Procedure. 

To check our variational approach, we have carried out 
an independent minimization to calculate the. ground- 
state energy of the system containing tt bondsEa, without 
making any assumptions about the presence or absence 
of fractional vortices. To carry out this minimization, we 
imagine that the ij bond is actually an overdamped 
Josephson junction connecting nodes i and j. The cur- 
rent flowing through that bond from node i to node j 



FIG. 10. Total energy as function of A for a single n bond 
obtained using (a) variational ansatz (*), and (b) numerical 
simulations (solid line). 

is then 

hi = Ic,ij sin {(pi -<t>j) + Yr~~^ ( d i ~ e i) j ( 37 ) 

and the sum of these currents must equal the total ex- 
ternal current, If xt , fed into node i: 

Y J kj = If xt - (38) 

b 

where J Ci y is the critical current of the junction between 
grains i and j, and Rij is the corresponding shunt re- 
sistance. These equations can be put into dimensionless 
form using the definitions iij = Itj / I c and gij = R/Rij, 
where I c and R are a convenient normalizing critical cur- 
rent and shunt resistance, and introducing the natural 
time unit r = H/ (2eRI c ). Combining these equations 
yields a set of coupled ordinary differential equations 
which is easily reduced to matrix form and solved nur 
merically, as described by many previous investigators^. 
For our work, we employed a fourth-fifth order Runge- 
Kutta integration with variable time step. 

For present purposes, we are interested, not in examin- 
ing the dynamical properties of arrays with tt bonds, but 
rather in finding the minimum-energy configuration of 
such arrays. To that end, we have simply iterated this set 
of coupled equations of motion, with no external current, 
allowing the phases to evolve until they reach a time- 
independent configuration. As has been shown by previ- 
ous workers, this configuration will correspond to a local 
minimum-energy state of the corresponding Hamiltonian 
H = — (hlc;ij '/2e) cos(0i — 4>j). We then compare 
the resulting configuration and energy with those pre- 
dicted by the fractional vortex variational ansatz for the 
ground state. 
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FIG. 11. Same as Fig. ^|but for two adjacent tt bonds. 



FIG. 12. Same as Fig. [ji] but for two adjacent tt bonds. 



To make the comparison as straightforward as possi- 
ble, we made the simplifying assumption iiy = R for all 
Josephson junctions, whether or tt. We took i Cj ij — 1 
for all normal junctions, and i c -ij = — A for all tt junc- 
tions. Since no external current is to be applied to the 
system, we carry out these calculations using square ar- 
rays of junctions with periodic boundary conditions in 
both directions. 

Each simulation begins with phases randomized at 
each gra in. The system is then relaxed according to the 
Eqs. ( ]37| ) and ( |38| ) for an interval of 50 — lOOr. We then 
evaluate the final energy, as well as the phase difference, 
8ij, across each tt junction. Once equilibrium is reached 
for a given A, we increment or decrement A, and the sys- 
tem is allowed to relax again without re-randomizing the 
phases. Even quite large arrays (50 x 50 plaquettes) re- 
lax quite quickly using this procedure, except near the 
critical point, but care must be taken to avoid taking 
data from simulations in which the system is trapped in 
a metastable state. We have used arrays ranging from 10 
x 10 plaquettes to 50 x 50, and occasionally as large as 
70 x 70 to examine convergence of equilibrium values. 

Figs. [| and |l^ show the the exact ground-state en- 
ergy and corresponding q for the case of a single tt bond 
in a host of normal bonds, as calculated by this numer- 
ical procedure. The results are also compared to the 
total energy obtained from a ground-state configuration 
corresponding to a pair of bound fractional vortices of 
charge ±g, calculated numerically. As shown in the fig- 
ures, the agreement is excellent, thereby indicating that 
the ground-state energy is indeed well characterized by a 
bound pair of fractional vortices. 

Figs, [n] and |l2| show a similar comparison for the case 
of two tt bonds. Once again, the results obtained numer- 
ically from the RSJ equations for both the total energy 



and the bond angle across the tt bonds, are in excellent 
agreement with those found from the fractional vortex 
ansatz, suggesting that the ground state, in the case of 
two 7r bonds, is again well characterized, over a range of 
A, by a pair of oppositely charged fractional vortices. 

Finally, we briefly discuss the accuracy of our varia- 
tional approximation for the "wave function," i. e., the 
phase distribution in the ground state. As is well known, 
a variational wave function may give an excellent value 
for the ground state energy, but a less accurate picture of 
the ground state configuration. In particular, our varia- 
tional approach ignores the spin- wave degrees of freedom 
in characterising the ground state of the frustrated XY 
model, but it is possible that they may be required to 
get an accurate description of the phase distribution. To 
test the accuracy of our variational phase distribution, we 
have compared it to the exact (numerical) phase distribu- 
tion in the ground state in several cases. The difference 
between the two configurations is shown graphically in 
Fig. 13 for the case of two tt bonds. As can be seen, the 



difference between the variational and exact wave func- 
tions is almost always less than 2-3% of the bond angle at 
the tt junction. We have looked at the results for one and 
two tt bonds for varying bond strengths, and in all cases 
considered the discrepancy is small. Thus, we conclude 
that the phase distribution as well as the energy is well 
approximated by our variational ground state involving 
only fractional vortices. 



IV. SUMMARY AND POSSIBLE SIGNIFICANCE 

The original motivation for this work was to study,* 
bonds in relation to the experiments of Kirtley et at^tB 
on ir-grain boundaries in high-T c superconductors. 
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FIG. 13. Graphical representation of the difference be- 
tween the variational phase configuration (two oppositely 
charged fractional vortices) and the numerically obtained 
ground state phase configuration for the case: A = 10, N n — 2 
(string of two adjacent parallel tt bonds). The phase differ- 
ence, as a fraction of the bond angle across the tt junctions, is 
shown on a gray scale given at the right edge of the diagram. 
A 70x70 lattice is considered and alternate stripes along the 
y-axis represent vertical and horizontal bonds. The tt bond 
string is shown in center of the Figure. 

Sigrist and Ricel have shown that the Josephson coupling 
across a grain boundary between two d-wave supercon- 
ductors can have either sign, depending on their crys- 
tallographic orientations, thereby giving rise to the pos- 
sibility of 7r-grain boundaries. Of course, in the present 
model, we are treating not a 7r-grain boundary, but rather 
a string of tt bonds in the discrete XY model. Never- 
theless, we argue that this string could be viewed as a 
crude model of such a grain boundary. It has been argued 
that even so-called 'single-crystal' high-T c superconduc- 
tors can be effectively represented as an array of super- 
conducting grains weakly interacting via the Josephson 
coupling between themc3. The typical lattice spacing for 
the high-T c materials in such a model has been quoted 
to be as large as 1 /im. Thus, the chain of tt bonds in our 
model can be taken as representing the coupling of grains 
across a 7r-grain boundary, and the length of the chain 
will depend on the dimensions of the grain boundary, and 
the interpretation of the effective lattice spacing. 

Now we turn to a summary of the experiments. The 
relevant experiments fall into two categories. In the 
tricrystal experiments, the intersections of three: grain 
boundaries at a "tricrystal point" were studiecollj. At 
special orientations of the grain boundaries, these ex- 
periments found that a half quantum of flux is trapped 
around the tricrystal point - a result that has been in- 
terpreted as verifying the cZ-wave symmetry of the super- 
conducting order parameter. In the tricrystal geometry, 
observing a trapped half-flux quantum can then be ex- 



plained by the fact that one ofj-the grain boundaries can 
be taken to be a 7r-boundaryBO. In the second class of 
experiments, a triangular (or a hexagonal) single-crystal 
superconductor was inserted into a single crystal super- 
conducting host of the same material, but with crystal 
axes misoriented with respect to those of the inclusion. In 
these systems, Kirtley et a£3 have found evidence of frac- 
tional (not half-integepl flux entrapment. These results 
have been interpreted! 2 ^ as evidence that the supercon- 
ducting order parameter violates time-reversal symmetry, 
either in the bulk or at an interface. Indeed, recent ex- 
periments have reported fractional fl-upc entrapment even 
in the absence of 7T-grain boundaries^!, possibly support- 
ing the existence of an order parameter which violates 
time-reversal symmetry. 

If, in the triangular inclusion, only one of the three 
boundaries is a 7r-boundary, the two "zero" boundaries 
will have little effect on the arrangement of the order 
parameter phases, and can reasonably be ignored. Sim- 
ilarly, in the tricrystal, if only one of the three grain 
boundaries is a 7r-boundary, this boundary would cor- 
respond to a semi-infinite chain of tt bonds, while the 
other two "zero" boundaries can again be ignored in the 
model. Thus, a finite chain of tt bonds may be suitable 
for modeling the triangular inclusions, and the extrapo- 
lation for long chain lengths is relevant for the tricrystal 
experiments. 

Next, we speculate about the relationship of our re- 
sults to the observed trapping of non-half-integers of flux 
in triangular inclusions. The trapped flux is usually re- 
lated to the phase differenjee across the grain boundary 
by the following argument!], which we restate to apply to 
our geometry. Consider a closed integration contour C 
(of radius r 3> a) centered at one end of the grain bound- 
ary, and passing through the grain boundary. We wish 
to consider the flux enclosed by this loop. The path is 
taken to be deep inside the grains, so that the Meissner 
effect dictates that the supercurrent density j = 0. Since 
j oc V(j> — (27r/<I>o)A, where <j> is the phase of the super- 
conducting wave function, <f>o — hc/(2e) is the supercon- 
ducting flux quantum, and A is the vector potential, it 
follows that 

27T 

W - (39) 

Now let C = C 1 + C 2 , where C 1 is the part of the 
contour not including the grain boundary. In the ap- 
proximation that C 2 can be taken to be infinitesimally 
short, the integral J c2 A • d£ ~ 0. In addition, we have 
J c , 2 V0 • d£ = A(j>, the phase discontinuity across the 
grain boundary. But also <f> must be continuous around 
C, modulo 2tt. Combining all these conditions with Eq. 
@, we find that 

2tt 

A<t> = 2ttu $, (40) 

where n is an integer and $ is the flux enclosed by the 
entire contour. Thus, the flux enclosed by C is related 
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site index, i site index, i 

FIG. 14. The calculated phase jump, A, which corresponds to the flux measured by a SQUID of size 10 x 10 lattice spacings, 
for a chain of n bonds of length n and bond strength A. The x coordinate denotes the position of the center of the SQUID 
relative to the leftmost n bond. A is defined to be the sum, taken in the counterclockwise direction, of the phase discontinuities 
across those bonds intersected by the perimeter of the SQUID, which are either n bonds in the grain boundary, or lie along 
the extension of the grain boundary along the ±x direction, (a) n = 20, A = 0.25. (b) n = 40, A = 0.25 (full curve) and 1.0 
(dashed curve). 



to the phase defect across the tt junction in the loop. In 
particular, if A<p is a non-half-integer fraction of 2tt in 
the ground state, then the flux enclosed will also corre- 
spondingly be fractional. Hence, a non-half-integer frac- 
tional flux is correlated with a phase jump across the 
grain boundary which is a non-integer fraction of tt. 

Our results show that A<p ^ tt for an interior bond 
in a finite chain of tt bonds. In fact A<j> can be 'tuned' 
to be any fraction of tt by simply varying the strength 
of the bonds for any finite chain length. Thus, a nec- 
essary condition for the occurrence of a non-half-integer 
flux quantum is indeed satisfied. But this result still 
does not demonstrate that the the trapped vortices cor- 
respond to non-half-integer flux quanta, because our cal- 
culations do not include the magnetic fields induced by 
the currents near the 7r-grain boundary. These fields will 
change A, and hence, the phase arrangement itself to 
some extent. Thus, we cannot rigorously infer the mag- 
netic flux when|-thcsc inductive effects are omitted from 
the calculations^. 

Although our present calculations do not include these 
inductive effects, it is still instructive to look at the phase 
distribution as if Eq. ( 40 ) were valid anyway. In particu- 
lar, let us try to model the flux configuration obtained by 
scanning the 7r-grain boundary using an idealized square 
SQUID. We take the flux through the SQUID to be the 
same as that through the corresponding contour C as de- 
scribed above. According to the argument just given, the 
flux passing through the SQUID is therefore proportional 
to the sum of the phase jumps A = A(j)i around the 
SQUID contour, across those bonds for which the phase 
has a discontinuity. In our simplified model for the flux 



through the SQUID, these discontinuities occur across 
the two bonds where the contour, taken counterclockwise 
around the SQUID, intersects the 7r-grain boundary or its 
extension along the x axis. 

In order to make a reasonable connection to the ex- 
perimental geometry, we estimate the lattice spacing a 
in our model using 



E.j = I c $o/c = a 2 J c $ /c, 



(41) 



where Ej is the Josephson coupling energy between ad- 
jacent grains, I c is the associated intergranular critical 
current, J c is the macroscopic critical current density set 
by the Josephson effect coupling, and a is the lattice spac- 
ing of the granular array. 

Using the experimental estimates for Ej and J c (see 
Ref. 13), we estimate a = l.l^tm, which iSpip agreement 
with the typical value for these materialseS. Since the 
triangular insertions are roughly 20/im in length (for each 
side), we have looked at the results for a tt bond chain 
of length 20 lattice spacings, with the SQUID diameter 
also taken to correspond to the experimental diameter. 
In Fig. 14(a), we show A as calculated for a tt bond chain 
of length n = 20 and strength A = 0.25, and a SQUID 
of diameter d = 10. In Fig. 0(b), we show the same 
for a chain of length n = 40. We note that for a fixed 
bond strength, A becomes more concentrated near the 
chain ends with increasing chain length n. Furthermore, 
we can also make A more localized near the chain ends 
by increasing A, as shown in Fig. fUjfb). 

The profile of A, shown in Fig. |l4|(a), strikingly re- 
sembles the flux profile measured in Ref. 13 across one 
side of a triangular insertion. In view of this similarity, 
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it would be interesting to see if the changes we see with 
bond strength and chain length are also found experi- 
mentally. (Experimentally, the strength of the bonds can 
be changed by varying the misoricntation angles of the 
inclusions.) If inductive effects do not change the qual- 
itative picture presented above, then these experiments 
would provide evidence supporting the interpretation of 
the grain boundary as a string of it bonds. .— . 

We also comment on the fact that Kirtley et a£3 were 
able to reproduce their measured flux configuration with 
a certain arrangement of fractional magnetic charges. 
Our picture suggests one way of understanding why this 
modeling works. The key is that the flux distribution is 
closely related to the gauge-invariant phase jump (Ajij) 
across the boundary, given by 

A 7ij = A^jj - (2i/$„)AA, J , (42) 

where i,j label sites connected by a bond across the grain 
boundary, and AAy and A<pij are the corresponding dis- 
continuities in the vector potential and the phase across 
that bond. A nonzero Ajij can therefore be attributed 
either to a nonzero A&j or a nonzero AAij (or a combi- 
nation of both) . In our calculations, we have assumed a 
nonzero Afaj and take AAij — across the branch cut. 
The opposite choice (AAij 0; Afaj — 0) corresponds 
to a fractional magnetic monopole, since the vector po- 
tential of a magnetiG-^nonopole changes discontinuously 
across a branch cutEl Since the physical quantity is 
the gauge-invariant phase difference, these pictures are 
equivalent. 

In summary, we have introduced a set of "fractional 
vortex " excitations in the XY model, and have derived 
an expression for the interaction energy of a bound pair 
of fractional vortices. Furthermore, we have studied the 
ground state of the XY model on a two-dimensional 
lattice containing tt bonds. For strings of n bonds of 
any length, we find that there exists a minimum bond 
strength, above which the ground state can be character- 
ized by pair(s) of oppositely charged fractional vortices. 
We have verified this ansatz by carrying out indepen- 
dent numerical simulations for the ground-state config- 
uration of this system. Finally, we have discussed the 
possible connection between these calculations and the 
trapped fractional flux quanta, which are observed near 
grain boundaries in high-T c superconductors. 
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